Numerical models of irrotational binary neutron stars in general relativity 
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We report on general relativistic calculations of quasiequilibrium configurations of binary neutron 
stars in circular orbits with zero vorticity. These configurations are expected to represent realistic 
situations as opposed to corotating configurations. The Einstein equations are solved under the 
assumption of a conformally flat spatial 3- metric (Wilson-Mathews approximation). The velocity 
field inside the stars is computed by solving an elliptical equation for the velocity scalar potential. 
Results are presented for sequences of constant baryon number (evolutionary sequences). Although 
the central density decreases much less with the binary separation than in the corotating case, it 
still decreases. Thus, no tendency is found for the stars to individually collapse to black hole prior 
to merger. 



OO 

On 



0< 



PACS number(s): 04.30.Db, 04.25.Dm, 04.40.Dg, 97.60. Jd 



o 

| Inspiraling neutron star binaries are expected to be among the strongest sources of gravitational radiation that 
could be detected by the interferometric detectors currently under construction (GEO600, LIGO, TAMA and Virgo). 
These binary systems are therefore subject to numerous theoretical studies. Among them are fully relativistic hydro- 
dynamical treatments, pioneered by the works of Oohara and Nakamura (see e.g. and Wilson et al. The most 
recent numerical calculations, those of Baumgarte et al. ||,[5j and Marronetti et al. ||, rely on the approximations of 
^vq • (i) a quasiequilibrium state and (ii) of synchronized binaries. Whereas the first approximation is well justified before 
| the innermost stable orbit, the second one does not correspond to physical situations, since it has been shown that the 
■ gravitational-radiation driven evolution is too rapid for the viscous forces to synchronize the spin of each neutron star 
with the orbit |?],[| as they do for ordinary stellar binaries. Rather, the viscosity is negligible and the fluid velocity 
circulation (with respect to some inertial frame) is conserved in these systems. Provided that the initial spins are not 
in the millisecond regime, this means that close binary configurations are well approximated by zero vorticity (i.e. 
irrotational) states. 

Moreover, dynamical calculations by Wilson et al. ^,|) indicate that the neutron stars may individually collapse 
into a black hole prior to merger. This unexpected result has been called into question by a number of authors (see 
£jrj! Ref. H for a summary of all the criticisms and their answers). As argued by Mathews et al. |J, one way to settle 
• • ■ this crucial point is to perform computations of relativistic irrotational configurations. We present here the first 
. ' quasiequilibrium irrotational relativistic binary neutron stars models computations. 

We have proposed a relativistic formulation for quasiequilibrium irrotational binaries ||To|| as a generalization of the 
Newtonian formulation presented in Ref. [TlJ| . The method was based on one aspect of irrotational motion, namely 
the counter-rotation (as measured in the co-orbiting frame) of the fluid with respect to the orbital motion (see also 
Ref. ) . Since then Teukolsky [jl3| and Shibata (l4|] gave two formulations based on the definition of irrotationality, 
which implies that the specific enthalpy times the fluid 4-velocity is the gradient of some scalar field [^5| (potential 
flow). The three formulations are equivalent; however the one given by Teukolsky and by Shibata greatly simplifies 
the problem. 

The hydrodynamic equations may be derived as follows. For a perfect fluid at zero temperature, the momentum- 
energy conservation equation V ■ T = is equivalent to the uniformly canonical equation of motion MlPffl , 



and 



u-(VAw)=0 (1) 



V ■ (nu) = , (2) 



where u is the fluid 4-velocity and w the momentum density 1-form w = hu, h being the fluid specific enthalpy 
h = (e +p)/(tob«) (V A w denotes the exterior derivative of w). In the above equations, n, e and p denote 
respectively the fluid proper baryon density, proper total energy density and pressure. It is clear that a potential flow 

w = V* (3) 
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is a solution of the equation of motion ([!]). Moreover, this particular solution is the relativistic generalization of the 
classical irrotational flow and, as stated above, corresponds to the physical situation reached by a binary system of 
neutron stars. 

As a first approximation of the relativistic treatment of the problem, we will assume that there exists a helicoidal 
symmetry p!o| . Let us denote by 1 the associated Killing vector. It is to be noticed that this symmetry is exact at the 
Newtonian limit. The helicoidal symmetry implies £[W = 0. From Cartan's identity Aw = 1 • V A w + V(l • w), the 
potential form ^ leads immediately to the following first integral of motion 

1 w = Const . (4) 

This was first pointed out by Carter jl7) . Note that this result is not merely the relativistic generalization of the 
Bernoulli theorem which states that 1 • w is constant along each single streamline and which results directly from 
the existence of a Killing vector without any hypothesis on the flow. In order for the constant to be uniform over 
the streamlines (i.e. to be a constant over spacetime), as in Eq. (|J), some additional property of the flow must be 
required. One well known possibility is rigidity (i.e. u colinear to 1) fil| . The alternative property with which we are 
concerned here is irrotationality. 

The fluid motion is now completely determined by the scalar potential ty. The equation for ^ can be derived from 
the baryon number conservation (^J). One obtains: 

^V- W + W- V (^) =0 . (5) 

h \hJ 

Within the 3+1 formalism and taking into account the helicoidal symmetry, this last equation is nothing but a 
Poisson-like equation which readsj^] 

nDiD 1 ^ + = — -^BWin + n < f + -^-BM A bxh - D^A In N - —Di(hT n ) I + KnhT n , (6) 

where we have introduced the covariant derivative D with respect to the spatial 3-metric, the trace K of the extrinsic 
curvature tensor, the shift vector B and the lapse N defined by the 3+1 orthogonal decomposition of the helicoidal 
Killing vector 1 = Nn — B (n being the unit future directed normal vector to the hypersurface t — Const) and the 
Lorentz factor r n = — n • u of the fluid with respect to the Eulerian observer whose 4-velocity is n. This equation has 
been recently derived by Teukolsky [|l3| and independently by Shibata |L4|]. Note that Eq. (|^) is independent of the 
gravitational field equations. 

As a first milestone in our project of studying coalescing binary systems, we will adopt the Wilson-Mathews 
approximation for the form of the metric |l^,§|. This approximation consists in taking a conformally flat 3-metric, so 
that the full spacetime metric reads 

ds 2 — ~(N 2 — B.B^dt 2 - 2Bidt dx i + A 2 r\ n ,dx l dx j , (7) 



where rj is the flat space metric. The field equations reduce now to the Wilson-Mathews equations [|2j|3j for TV, B 
and A. Right now, it is not obvious whether the Wilson-Mathews approximation is valid in the case of coalescing 
compact binaries. We presently use this particular form of the metric in order to simplify the problem. However, 
it is to be noticed that (i) the 1-PN approximation to Einstein equations fits this form, (ii) it is exact for arbitrary 
relativistic spherical configurations and (iii) it is very accurate for axisymmetric rotating neutron stars |20[ . An 
interesting discussion about some justifications of the Wilson-Mathews approximation may be found in Finally 
we chose maximal slicing: K = 0. 

Following (k]], we introduce Q such that the helicoidal Killing vector 1 satisfies 

where r and <p are respectively the time and azimuthal coordinate associated with the asymptotic inertial observer 
at rest with respect to the binary system (i.e. such that the ADM 3- momentum vanishes on the slices r = const). 
Besides, we introduce the non-rotating shift vector N defined by 



1 Latin indices run in the range 1,2,3 and geometrized units (G = 1 and c = 1) are used. 
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(9) 



The gravitational field equations are derived within the 3+1 formalism from the Hamiltonian constraint, momentum 
constraint and trace of the spatial part of the Einstein equation (|J^]. Introducing v = lniV and (3 — ln(AZV), they 
can be written as 

A/3 = AkA 2 S + ^A 2 K i0 K %1 - i (v^vV + V^V>) (10) 
Ay = An A 2 {E + S) + ,r-'/\,,A'' - WiiAT/3 (11) 

AN 1 + -V (VjN j ) = -16ttNA 2 (E + p)U l + 2NA 2 K ij W j (3/3 - Av) , (12) 
3 

where V is the covariant derivative associated with the flat 3-metric r) and A = V*Vi is the corresponding Laplacian^. 
E = r^(e + p) — p, S = 3p+ (E + p)UiU l and U l — D 1 ^ /(hT n ) are respectively the fluid energy density, trace of the 
stress tensor and fluid 3- velocity, the three of them measured by the Eulerian observer. T n can be computed according 
to 



r n 



1 + W V ^* 



1/2 1 



Vi - UiU i 

and the extrinsic curvature tensor is computed by means of the identity 



K ij = -— ^ {VaP + V j N< - ^V k N k J , 



(13) 



(14) 



which results from the Killing equation for 1. 

The matter distribution is determined by the first integral of motion (Q). Taking its logarithm leads to 



1, A , 2 B*B' 



H + v+-\n[l- A%,—) + In T = Const , (15) 

where H := ln/i and T is the fluid Lorentz factor with respect to the co-orbiting observer (i.e. the observer whose 
4- velocity is collinear to 1): 

r = r n - N . . . (i6) 

„ B I B J 

Note that Eq. ( p"5|) corresponds to Eq. (66) in Ref. and that the constant which appears at its r.h.s. is nothing 
but the logarithm of the constant C introduced by Teukolsky jlj| . 

Now, introducing Q = d In H/d Inn, the equation for the fluid velocity potential (^) becomes 
(HMr + THViV = -A 2 hT n ^ t H + (H j (V* + A 2 hT n ^-j - vW l( 3 - A 3 ^i(hT a )} . (17) 

The equations to be solved (0) (0) © (0) constitute a system of non-linear Poisson-like equations. Because 
of the elliptical nature of these equations, we have exhibited the common flat Laplacian operator A which can be 
solved by means of usual spectral methods (cf. e.g. 01 or [p2|). Because of the non-linearities, we use an iterative 
procedure based on a multi-domain spectral method Ba] to get the solution. The numerical code will be described 



2 Throughout the article, we use the notation V = rf^Vj 
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in details in a forthcoming paper j24|. Let us simply mention here some tests passed by the code. In the Newtonian 
and incompressible limit, the analytical solution constituted by a Roche ellipsoid is recovered with a relative accuracy 
of ~ 10~ 9 (cf. Fig. 6 of Ref. p3|). For compressible and irrotational Newtonian binaries, no analytical solution 
is available, but the virial theorem can be used to get an estimation of the numerical error: we found that the 
virial theorem is satisfied with a relative accuracy of 1CT 7 . A detailed comparison with the irrotational Newtonian 
configurations recently computed by Uryu & Eriguchi |25|]26| ] will be presented in j23] . Regarding the relativistic case, 
we have checked our procedure of resolution of the gravitational field equations by comparison with the results of 
Baumgarte ct al. Q] which deal with corotating binaries [our code can compute corotating configurations by simply 
setting lnT = in Eq. (|l^) and using U l — —B l /N for the fluid 3- velocity]. We have performed the comparison with 
the configuration za = 0.20 in Table V of Ref. ||. We used the same equation of state (EOS) (polytrope with 7 = 2), 
same separation rc and same value of the maximum density parameter g max . We found a relative discrepancy of 1.1% 
on Q, 1.4% on M , 1.1% on M, 2.3% on J, 0.8% on z A , 0.4% on r A and 0.07% on r B (using the notations of Ref. [§). 
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FIG. 1. Relative variation of the central energy density e c with respect to its value at infinite separation ej. nf as a function 
of the coordinate separation d (or of the orbital frequency Q./(2tt)) for constant baryon mass Mb = 1.625 Mq sequences. The 
solid (resp. dashed) line corresponds to a irrotational (resp. corotating) sequence of coalescing neutron star binaries. Note that 
there is no substantial increase of the central density as the separation decreases. 



These tests being passed, we turned towards calculations of irrotational relativistic binaries. We chose to investigate 
the instability issue raised by Wilson and Mathews Q by computing an evolutionary sequence (i.e. a sequence at fixed 
baryon number) made of irrotational quasiequilibrium models. We took the same configuration than that presented by 
Mathews, Marronetti and Wilson (Sect. IV A of Ref ||), namely two identical stars obeying a 7 = 2 polytropic EOS 
[p = K(mBn) 7 , e = m-Bii+p/ (7— 1)] with k — 1.8x 10~ 2 J m 3 kg~ 2 and each having a baryon mass Mb = 1.625 Mq. For 
such parameters, we found that the gravitational mass of a single star in isolation is M = 1.515 Mq (in agreement with 
Ref. H), with a central energy density e" lf = 4.005 p nuc c 2 (p n uc '■— 1-66 x 10 17 kgm -3 ) and a compactification ratio 
M/R = 0.140 (e c and M/R are slightly different from that quoted in Ref. || probably due to different values for the 
constants G, c, Mq and m B ; we use G = 6.6726x 1CT 11 m 3 kg~V 2 , c = 2.99792458 x 10 8 ms" 1 , Mq = 1.989xl0 30 kg 
and m B = 1.66 x 10" 27 kg). 

We define the coordinate separation d as the coordinate distance between the two density maxima. Using the same 
value as the one considered by Mathews et al. ||, namely d = 100 km, we found that a Mb = 1.625 M© irrotational 
configuration in quasiequilibrium at this separation has a total angular momentum of J / (2Mb) 2 = 1-13 which is quite 
similar to the value 1.09 found by Mathews et al. J^]. However, we did not observed any substantial increase of the 
central (i.e. maximum) density with respect to static stars in isolation, as they did (they report a central density 
increase of 14%). 

In order to investigate the evolution of a coalescing binary system, we have computed a full sequence at fixed 
M B = 1.625 M , starting at the separation d = 110 km (fi/(27r) = 82 Hz) and ending at d = 41km (fi/(2vr) = 332 Hz). 
We have considered both corotating and irrotational cases. The evolution of the central density along these sequences 
is shown in Fig. [TJ. In the corotating case, we found that the central density decreases substantially when the stars 
approach each other, as expected from previous independent calculations |^,^|. In the irrotational case, we found that 
the central density still decreases with the separation but much less than in the corotating case. 
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All the calculations have been performed with spectral methods using the library LORENE (Langage Objet 
POUR LA RElativite NumeriquE [p7j). We warmly thank Jean-Pierre Lasota for his support. The calculations 
presented in this paper have been made possible thanks to a special grant from the SPM and SDU departments of 
CNRS. 
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